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SUMMARY 


ANALYSIS OF RANDOM STRUCTURE-ACOUSTIC INTERACTION PROBLEMS 
USING COUP LED BOUNDARY E LEM ENT AND FINITE E LEM E N T METHODS 

A coupled boundary element (BEM) -finite element (FEM) approach 
is presented to accurately model structure— acoustic interaction 
systems. The boundary element method is first applied to interior, 
two and three-dimensional acoustic domains with complex geometry 
configurations. Boundary element results are very accurate when 
compared with limited exact solutions. 

Structure-interaction problems are then analyzed with the cou- 
pled FEM-BEM method, where the finite element method models the 
structure and the boundary element method models the interior acous- 
tic domain. The coupled analysis is compared with exact and ex- 
perimental results for a simplistic model. Composite panels are 
analyzed and compared with isotropic results. The coupled method 
is then extended for random excitation. Random excitation results 
are compared with uncoupled results for isotropic and composite 


panels . 


FORWARD 


This technical report contains the research results on random 
structure-acoustic interaction problems using coupled boundary 
element and finite element methods. The research is sponsored 
by a NASA Langley Research Center Grant NAG-1-1358, entitled 
"Nonlinear Analysis of Composite Aerospace Structures in Sonic 
Fatigue." Technical monitor is Dr. Stephen A. Rizzi, Structural 
Acoustics Branch, Fluid Mechanics and Acoustics Division, NASA 
Langley Research Center. The authors are extremely thankful 
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Chapter 1 


INTRODUCTION 

1.1 Preliminary Remarks 

Over the past few years, structural-acoustic interaction 
problems have shown enhanced interest in the aerospace and 
automotive industries. New supersonic aircraft designs such 
as the High-Speed Civil Transport (HSCT) , the National Aero- 
Space Plane (NASP) , and the Advanced Tactical Fighter ( ATF) 
have lead to some very complicated engineering problems. 
Acoustic excitation of the fuselage has been one of the main 
concerns during certain flight operations. In Figure 1.1, the 
acoustic loads on a B-52 wing during take-off are shown. The 
acoustic sound pressure levels reach as high as 164 dB [1]. 
Systems on new supersonic aircraft will be highly random in 
nature with temperatures possibly reaching 2000°F and acoustic 
sound pressure levels of 190 dB, as shown in Figure 1.2 [2] . 
In such extreme environments, the structural-acoustic system 
can be driven into an extreme random regime from such things 
as turbulent boundary layers. Under such conditions, the 
coupling of the acoustic and the structure systems and 
determining their combined random responses become extremely 
important . 
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Due to the harsh environments generated during flight, 
structural design becomes critical and strength-to-weight 
ratios also become vitally important. In extreme 
environments, composite materials can have major advantages 
over isotropic materials in the random response of structural 
components. Composite materials offer adjustable design 
variables such as lamination angles, layer thicknesses and 
fiber/matrix composition ratios. 

Due to the complex geometry of many structural problems, 
numerical methods have become the tool of choice. The finite 
element method (FEM) is a well established technique that 
offers many advantages when modelling the structure. The 
boundary element method (BEM) has become a popular technique 
when modelling acoustic domains. Recently, some attempts have 
been made to couple the boundary element method and the finite 
element method for solving structural-acoustic interaction 
problems. By coupling the two numerical techniques, the 
designer creates a very powerful tool for modelling real 
system behavior, including random response. 

1.2 Review of Previous Work 

Over the past 15 years, research on the indirect boundary 
element method (IBEM) and the direct boundary element method 
(DBEM) has increased. The DBEM solves directly for the 
acoustic pressures and velocities. The IBEM creates source 
potentials to model the boundary as acoustic sources and from 
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the potentials, acoustic pressures and velocities are 
obtained. The direct method is better suited to couple with 
the finite element method. The major differences between the 
IBEM and the DBEM have been documented [3-5]. 

Approximately, 80% of all boundary element papers dealing 
with acoustics are in the radiation and scattering area [6- 
14]. Acoustic radiation and scattering has also become an 
important research area for acoustic systems. The majority of 
research papers concerning internal domains have been in 
three-dimensional space. Two-dimensional problems present 
difficulties associated with the singularities of the Hankel 
function [15-19]. Research in three-dimensional duct systems 
has been conducted in three main areas. Room acoustics with 
acoustic sources started three-dimensional BEM research [20- 
22]. Another area of research has been in the design of 
mufflers. Here, the boundary element method has resulted in 
a very useful modelling tool for solving difficult geometry 
and boundary condition problems [23-26]. Research interest 
has also been directed towards the design of automobile 
cavities. The interior of automobiles can give rise to some 
very difficult geometry regions as well as reflecting and 
absorbing material boundary conditions [16,20,27,28]. 

Additional research in acoustics using the boundary 
element method includes element design and time domain 
solutions. Many types of boundary elements have been 
developed and used for acoustic domains. The first type of 
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elements used for acoustics were the constant, linear and 
quadratic elements [29], Triangular and quadrilateral 
elements have become increasingly popular with boundary 
element users [5,10,30,31]. Limited research has been 
conducted on transient problems with time domain solutions 
[32,33]. Coupling the finite element and boundary element 
methods has been a main focus of the research community over 
the past ten years. The initial idea to couple the 
discretized equations of FEM and the boundary integral method 
was developed by Zienkiewicz, Kelly and Bettess in 1977 [34]. 
The boundary element and finite element methods have been 
coupled for many types of problems. Numerous publications 
have linked the two methods in structural dynamics. Crack 
propagation [35], soil-plate coupling [36], and beam and plate 
vibration systems are some typical examples. Specific beam 
and plate problems include studies on structures with cut-outs 
[37,38], varying BEM and FEM regions [39-42], and circular 
plates with slots [43]. The coupled BEM-FEM technique has 
been used to model fluid-filled tanks and spheres used in the 
Space Shuttle [44,45]. Coupling FEM and BEM for the Laplace 
equation has shown limited interest [34]. 

Modelling structure-acoustic interaction problems has 
also been at the forefront of current research over the past 
few years. Radiation and scattering problems have received 
much of the research attention, since the boundary element 
method handles infinite regions very well [46-51]. Yuying 
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studied acoustic radiation of floating plates with odd 
geometries [52]. Coupling of interior acoustic and structure 
problems has also been popular. In 1984, Suzuki, Imai and 
Ishiyama provided the initial formulation coupling the modal 
FEM and BEM [53]. Many articles use a well-known cavity- 
backed plate problem to compare results [53-56] . The cavity- 
backed plate problem was first analyzed by Guy and 
Bhattacharya in 1973 [57]. They provided valuable semi- 
analytical results using modal approximations and also 
experimental data. In 1992, Bokil derived an exact solution 
to the coupled cavity-backed plate problem using Laplace 
transforms [58]. After the coupling technique showed that it 
was accurate and versatile, many interesting real world 
problems have been modelled. Aircraft fuselage studies were 
conducted in 1991 by Fyfe, Coyette, and van Vooren [59]. Many 
authors have attempted to model complex geometries of 
automobile cavities [59-61]. Some automobile cavities have 
used complex boundary conditions including fiber and foam 
absorbing material [20,60]. Other interesting automotive 
problems have been studied such as engine block noise [59] and 
gearbox noise [62] reduction studies. 

1.3 Objectives and Scope 

The overall objective of the present study is to present 
a coupled boundary element and finite element technique to 
solve coupled structural-acoustic interaction problems for 
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harmonic and random response. To perform this task, four 
steps are required. The first objective is to develop and 
validate the boundary element method as an accurate numerical 
technique for acoustic domains. The second step is to acquire 
a finite element program to accurately model the structural 
system. The third objective involves coupling the boundary 
element method and the finite element method to model the 
total coupled system. Coupling the acoustic and plate 
displacements as well as the acoustic force becomes very 
important. The first three steps have been completed 
previously. The new contribution of this study includes the 
coupled BEM-FEM technique to model realistic problems, 
including composite materials and random responses of the 
acoustic and structural domains. 

The boundary element method is well suited for interior 
and exterior acoustic problems. BEM reduces the 
dimensionality of the discretized domain by one dimension, 
thus making it very appealing for complex geometric domains. 
For instance, a three-dimensional domain is modelled by two- 
dimensional surface elements only. The two-dimensional 
acoustic problem is analyzed first. Two-dimensional domains 
are very difficult to model with BEM because of the 
singularities that occur in the assumed free— field pressure 
solution. Since the assumed pressure, or Green's function, 
contains the first and second order Bessel functions, the 
domain integral becomes singular at certain positions. This 
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study derives a new method of approximating these singular 
integrals. A few two-dimensional cases are modelled using 
linear and constant boundary elements including: mufflers, 
ducts with sudden area changes, and a three-section duct 
system. BEM results are compared with exact solutions for 
certain rectangular duct systems in Chapter 4. 

Three-dimensional duct systems are also modelled in this 
study. Singular integrals are easily evaluated by a change of 
coordinate systems. BEM results are again compared versus 
exact solutions for rectangular duct domains in Chapter 4. 
Irregular-shaped ducts are also modelled to show the 
versatility of the boundary element method. Chapter 3 
includes the coupled formulation for BEM and FEM. A brief 
introduction to the finite element formulation is given along 
with pertinent element-type information. Coupling techniques 
between BEM and FEM for random analysis is derived in Chapter 
3. Chapter 4 includes the case results from each system 
listed above. Results include two and three-dimensional duct 
systems and coupled structural-acoustic interaction problems. 
Random effects are also shown for the coupled interaction 
problem. Concluding remarks and recommendations for future 
work are presented in Chapter 5. 
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Chapter 2 


BOUNDARY ELEMENT FORMULATION 

In this chapter, the governing equations are derived for 
two and three-dimensional acoustic domain problems. The 
system equations and solution procedures are also given. 

The derivation begins with the wave equation and 
progresses with the derivation of the influence matrices, [H] 
and [G] . There is no stream flow considered and only the 
acoustic pressure is assumed present. Throughout this study, 
only linear acoustic and structure theory are used in the 
formulation. 


2.1 Acoustic Wave Equation 

The derivation of the linear wave equation is well 
defined and can be found in any acoustic text. The acoustic 
medium is assumed to be homogeneous, isotropic and perfectly 
elastic. The perfectly elastic assumption allows us to use 
the particle displacements and velocities of acoustic waves in 
the governing equations, similar to elastic waves in solids. 
No viscous or gravitational effects are considered in the 
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derivation. The equation of state of a perfect fluid medium 
is given as 


P = Bs, (2.1) 

where B is the adiabatic bulk modulus and s is the 
condensation. Only small displacements are considered for the 
acoustic particle displacement and the acoustic pressure, P, 
is considered to be small compared with the equilibrium 
pressure. The equation of continuity can now be applied to a 
volume element. The linearized continuity equation is shown 
as 


+ V- v a = 0, (2.2) 

where v, is the particle velocity and V is the divergence 
operator. 

The next step is to apply Newton's second law to the 
equilibrium of forces on the acoustic particles, and thus the 
inviscid linearized Euler's equation is produced 


Pa 


dV a 

!Tt 


-VP, 


(2.3) 


where p, is the density of the acoustic medium. Combining 
Equations (2.1), (2.2) and (2.3), and eliminating the 

condensation, s, yields the linearized wave equation 
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(2.4) 


v 2 p= _L 


d 2 p 

at 2 


where c is defined as, c=(B/p,) m . A detailed derivation of 
the linear wave equation can be found in Kinsler and Frey 
[65]. The pressure in Equation (2.4) is represented as 


P = P(x,y, t) . 


(2.5) 


2.2 Two-Dimensional Formulation 

Very little work has been done in the two-dimensional 
study of the boundary element method due to the difficulty in 
calculating the integrals of the fundamental solutions. V 2 , 
from Equation (2.4) , is the two-dimensional Laplacian operator 
given by 

v==JL t J!_. 

dx 2 by 2 

The fundamental acoustic pressure response is given as 

P=p(x,y)e tul (2.6) 

and thus the two-dimensional linear wave equation reduces to 
the Helmholtz equation, 

V 2 p+k 2 p= 0 , (2.7) 

where k is the acoustic wavenumber and is represented as 
k=Gj/c. The value o>, is the acoustic input frequency. The 
variable p represents the pressure distribution and q is the 
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partial derivative of the pressure with respect to the normal 
or, 

< 3 = -j^ = ~iPa (j}V a’ ( 2 . 8 ) 

Equation (2.8) is a direct result from Euler's equation as 
seen in Equation (2.3). This condition will be utilized as 
the boundary condition of the sound pressure on the surface of 
the duct, where d/dn is the outward normal derivative. 

Noting that p is only a function of spatial coordinates, 
the method of weighted residuals statement of the Helmholtz 
equation is given by 

Jjv 2 p+/c 2 p)p * dn =0 , (2.9) 

where the fundamental solution, p*, is the free-field Green's 
function for the Helmholtz equation and satisfies 

V 2 p* +k 2 p* +a=0, (2.10) 

in the domain n and a is the Dirac delta function. The 
fundamental solution for two-dimensional space and solution to 
Equation (2.10) is 

p-=±Hl l) (kr) , (2.11) 

and 
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(2 . 12 ) 


g* =— ( kr ) cosfl, 

where r is the distance between the source point, S, and the 
observation point, B, as seen in Figure 2.1 and H(kr) is the 
Hankel function [29]. The variable 0 is the angle between the 
normal of the j* element (observation point, B) and the 
distance vector, r, which is between the source point, S, and 
the observation point as shown in Figure 2.1. Figure 2.1a is 
a model of a two-dimensional duct utilizing constant boundary 
elements. Figure 2.1b is a two-dimensional duct utilizing 
linear boundary elements. Note that both ducts have inlet and 
outlet boundaries as shown. 

The Hankel function is defined as 

Hi u (kr) =J a (kr) +iY c (kr) , 


and 

H\ l) (kr) =J, (kr) +iY x (kr-) , (2.13) 

where J and Y are the Bessel functions of the first and second 
kind, respectively [66]. The subscripts on the Bessel 
functions can be identified as: 0 representing the zero order 

and 1 representing the first order. The definitions for the 
Bessel functions of the zero and first order are given as 

J-Yjx) = -Y^x) , jLjJx)=-J l (x) . (2.14) 
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Outlet 





is the Hankel function of the first kind and zero order 

and is the Hankel function of the first kind and first 

order . 

Integrating Equation (2.9) twice by parts and collecting 
terms, the governing equation becomes. 



-| r gp*dr+| r pg*dr, 


(2.15) 


where r represents the total boundary. The domain integral 
term in Equation (2.15) is evaluated by introducing the Dirac 
delta function from Equation (2.10), and thus reduces the 
integral to -c m p m . The c m term is constant for each boundary 
element and is defined as, 


c m =_ on the boundary (T) 

c =0 outside the boundary (n') ... 

( 2 * X 6 j 

c m =X on a smooth boundary (T) 
c m =l inside the boundary (fl) 

where a' is the internal angle, in radians, between connecting 
elements and the subscript m represents the source element. 
For example, a smooth element has a c m value of 1/2, since the 
angle at the representing node is tt. For smooth constant 
elements, the angle is always n , since the representing node 
is at the center of the element. For linear elements, a 
corner element has a tt/ 2 angle and a c m value of 1/4. For 
linear elements, a c m value is given for each node of the 
element. Figure 2.2 clearly shows the definition of the angle 
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a' for constant and linear elements. A complete derivation 
and description of the c m value is given by Brebbia [29]. By 
introducing the Dirac delta function into the domain integral, 
the domain integral is reduced to a constant term and thus the 
governing equation reduces to only boundary expressions. This 
assumption of using the free-space Green's function is the 
backbone idea of the boundary element method. 

Discretizing the boundary into elements, the governing 
equation can be written as, 

P m + E ( r <7*P dr= £ L P*<? dr, (2.17) 

where N is the total number of elements and j represents the 
j* element. The fundamental solutions from Equations (2.11) 
and (2.12) can now be substituted in Equation (2.17). At this 
point, one must decide which type of element to use. 
Extensive work has been completed using three types of 
elements: constant, linear, and quadratic elements. For the 
two-dimensional analysis, only constant and linear type 
elements are discussed. However, higher order elements can 
easily be derived and programmed. This study was not an 
attempt to employ higher order elements. Rather, the 
objective of this study is to understand and apply the 
governing equations of the boundary element method to model 
acoustic domains. 
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2.2.1 Constant Elements 


If constant elements are used, p and q are assumed 
constant over each element. The element is represented by a 
line and a node at the center of the element as seen in Figure 
2.3 below. 

Node Represents Element 

I 

© 

1 

Figure 2.3 Two-Dimensional Constant Element. 


Therefore, p and q can be extracted from the integral terms in 
Equation (2.17). Substituting p' and q*, from Equations (2.11) 
and (2.12), into Equation (2.17) and using constant elements. 


c p +TP 

m Sr m / - 
>1 


±±H\ l) (kr) 


N 


cos 9 dr p,=£ 


, 


±H?(kr) 


k?r 


g,.. (2.18) 


Note that Pj and qj have been pulled outside of the boundary 
integrals. Also, note that the angle 9 varies as the integral 
is enforced over the element. Equation (2.18) is the final 

form of the governing boundary element equation for constant 
elements. 
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2.2.2 Linear Elements 


For linear boundary elements, p and q are defined in 
terms of their end nodal values and two interpolation 
functions and $ 2 . The linear element is represented by a 
line with nodes at each end of the element as shown in Figure 
2.4. 

Both. Nodes 
Represent Element 

/ 

• 

1 2 
Figure 2.4 Two-Dimensional Linear Element. 



The interpolation functions are 

*' * 

and 

*2 = |(l + £), ( 2 . 19 ) 

where | is a dimensionless coordinate varying from -1 to +1. 
Now p and q can be expressed as 

pa) = *,p’ + *2 p 2 = [*i j^j (2,20) 


and 
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g(f) = *#' + * 2 g 2 = [«i *2] j^j- 


( 2 . 21 ) 


The values p 1 , p 2 , q 1 , and q 2 in Equations (2.20) and (2.21) are 
the constant pressures and normal derivatives of pressure at 
nodes 1 and 2, as shown in Figure 2.4. Substituting the 
expressions for p and q from Equations (2.20) and (2.21) and 
the fundamental solutions for p* and q* from Equations (2.11) 
and (2.12) into the governing boundary element equation, the 
governing equation for linear boundary elements becomes 


c p 

m c'm / v 

>1 


j r [*. * ! |-^H!' , (kr)]cos» drjj^J 


2N . 

-E j,[*. 


dr 


q[ 

g 1 

q\ 


( 2 . 22 ) 


The right side of Equation (2.22) is summed from 1 to 2N 
because q is the variable of the integral. In other words, 
there is one q value coming into the node and one q value 
going out of the node. So one boundary node has two q values, 
one q in and one q out, as shown in Figure 2.5. The subscript 
on q represents the node number and the superscript represents 
incoming or outgoing, (l=out, 2=in) . Accounting for the 
difference in q values, requires considerable more 
computational time and effort than constant elements. In 
Equation (2.22), note that only the #, and $ 2 values vary with 
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coordinate location and the p's and q's are constant nodal 


values . 


Boundary Node 



q 1 q 2 q 1 q 2 q 1 q 2 

1 1 2 2 3 3 

Figure 2.5 Definition of q Values for Linear Elements. 


2.3 Three-Dimensional Formulation 

Over the past few years, a vast amount of research has 
been performed on three-dimensional acoustic problems using 
the boundary element method. Most of the research has been in 
the acoustic scattering and radiation areas. More research 
has been done on three-dimensional domains than two- 
dimensional domains because two-dimensional problems involve 
Bessel functions and three-dimensional problems involve 
exponential functions. Thus integration solutions are much 
easier to work with in three-dimensions. 

The three-dimensional formulation is very similar to the 
two-dimensional formulation. The same governing boundary 
element equation exists for three-dimensions, 
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(2.23) 


C * Pm + E L V'P dr= E f p P*^ dT ' 

where the c m value is evaluated from Equation (2.16) in the 
same manner as the two-dimensional analysis. Since constant 
elements were found to give very accurate solutions compared 
with linear elements, only constant elements are used in the 
three-dimensional analysis. Therefore, p and q can be 
extracted from the integrals in Equation (2.23). A three- 
dimensional constant element is defined by four boundary nodes 
and represented by a node at the center of the element as seen 
in Figure 2.6. 


Xode Representing 
Element 



Boundary Node 


Figure 2.6 Three-Dimensional Constant Element. 
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The fundamental solution for three-dimensional space is 


and 


p* =e' ikr /4Tir 


q'=~— (l/?+i/c/r)cos0. (2.24) 

4tt 

Substituting Equation (2.24) into Equation (2.23), the 
governing three-dimensional boundary element equation becomes 


N 


Pm' E 


h i 


r N 

| fL— (l/P+i.k/r)cos0 dr p,=y^ |j e'^'fiirr drj g^, (2.25) 


where the total boundary has been discretized into N elements. 
The values Pj and qj have again been pulled outside the 
boundary integrals. Note that if linear elements were 
employed, the same approach as the two-dimensional case would 
follow with two dimensionless parameters, one in the x- 
direction and one in the y-direction. 


2.4 System Equations and Solution Procedures 

By applying the governing boundary element equation to 
each element, a set of equations is formed which can be 
represented in matrix form. The two-dimensional integrals in 
Equations (2.18) and (2.22) can be written as, 

K ] = ( r ;T tfl<1)( ^ )cosfl dr 


and 
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K ]=f r j H “’ ( * r)dr 


(2.26) 


for constant elements, and 



if ik „(i) 
2j — H \ 


(*r) 


cos 6 dr 


and 




<?r 


(2.27) 


for linear elements. The subscript m represents the source 
element and j represents the observation boundary element. 
The three-dimensional integrals in Equation (2.25) can be 
expressed as 


and 


[**] — 


-k 


— — (l/?+ik/r)cosfl dr, 
4 71 


(G m> ] = f e-“747rr dr. (2.28) 

J 1 / 

From this point, the solution procedure is the same for 
two and three-dimensions using the appropriate [H^] and [G mj ] 
matrices. Using matrix form, the governing boundary element 
equations can be written as 

*.p. } -[«-]{»«}• < 2 - 29 > 

The subscript m sums from 1 to N, j sums from 1 to N and k 
sums from l to N for constant elements and from 1 to 2N for 
linear boundary elements. 
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The c m p m term can be combined with the [fi^] matrix to 
form a new matrix, [H^] , where [H^] is given as 


H mi = 

mj 


H • + C 
H 

mj 


m=j 

m?£j 


( 2 . 30 ) 


Therefore, the system of equations in matrix form is, 


[ H ] { p } = [G] {g}. 

The matrices [H] and [G] are called the influence matrices 
because they relate how one element influences another. The 
integral expressions in the influence matrices are solved with 
a four-point Gauss integration scheme as shown in Appendix A. 
Both {p} and {q} vectors contain unknown values of p and q, 
respectively. For a particular element, p and/or q may be 
unknown. If an impedance boundary condition exists, then p 
and q are unknown and only the relationship between p and q is 
known. To solve Equation (2.31), the matrices must be 
transformed into a new matrix equation 

[A] {X} = {B}, (2.32) 

where all of the unknown p and q values are collected in the 
vector {X}. The matrix [A] and the vector {B} contain complex 
values and hence all unknown p and q terms can be solved on 
the boundary. Also note that the matrix [A] is fully 
populated and unsymmetric. If any pressures or velocities of 
internal points are required, they can be found as a post 
processing operation only after the boundary pressures and 
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velocities have been calculated. The same system equation is 
used for determining pressure and velocity with the 
coordinates of the internal point and a different c m value, 
substituted into Equation (2.29). 


2.5 Boundary Element Singularities 

Evaluating singularities is a big concern when using the 
boundary element method. Singularities occur when the 
observation point lies on the source element. Usually, the 
term does not produce any complications for smooth 
elements, since the normal is perpendicular to the element. 
When this occurs, H,^ equals zero due to the cos0 term. 
However, the G,,,,,, term almost always yields a singularity. 


2.5.1 Two-Dimensional Problems 

For the two-dimensional Helmholtz equation, a singularity 
occurs in the integral for constant elements, where the G,,,,,, 
integral term is, 


G ™" = Jr ^ H o^ kr ) dr = ^\ r { J o(kr) + iYJkr))dr. (2.33) 

This study provides approximations to the singular integrals 
of the Bessel functions. As r approaches 0, the Y 0 (kr) term 
tends to -<®. The integral in Equation (2.33) can be broken 
into two independent integrals, 


Gmm 4k jo J °( e J de ° + 4/fjo Y oi € o) d€ o‘ 


(2.34) 
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A variable change was performed to help with the integration. 
The first integral is given by, 


^J 0 (e 0 )de 0 =klJJkl) 

+ ?*L[s n (kl)J 1 (kl) -S T 1 (kl)J 0 (kl)], (2.35) 

where S T0 and S T1 are the Struve Functions and 1 is the length 
of the element. The second integral is expressed as. 


( e a ) de 0 =klY a (kl ) + [s ro (kl ) Y, (kl ) -S n (kl ) Y a (kl ) ] , (2.36) 

where the exact integration above is given in [66]. 

For linear elements, the G,^ term also produces a 
singularity. The integral term is, 

= J p *i H< o U (*?> dr = j j r ( i-O [ J A k r) +iY 0 (kr) ]dr . (2.37) 

By changing the coordinate system from r to /3 0 , the initial 
integral is broken into four minor integrals, 


G = 

nun 


— f“ 

8 k Jo 


28 28 i. 

2 JM + 2 1Y O (0,) --£j 0 (( 3 a ) --rr °—Y 0 (& 0 ) 


kl 


kl 


W n 


(2.38) 


The first two integrals in Equation (2.38) are solved using 
Equations (2.35) and (2.36), respectively. The third integral 
is solved with the identity, 


PJM Wo = kU 1 (.kl) , (2.39) 

and the fourth integral is solved by, 


v 
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J ki 2 r ( i ) 

0 o Y o (0 o )d0 o = klY t (kl)* l , (2.40) 

0 7T 

where r f is the Gamma function. A small error is produced in 
the terms, because the Struve functions are calculated by 
a series expansion approximation. The complete definitions 
and approximations of the Struve functions are given in 
Appendix A. Since the two-dimensional singularities are such 
a problem, many authors did not address two-dimensional 
acoustic problems. Working through the singularities is a 
difficult obstacle to overcome when using the boundary element 
method . 

2.5.2 Three-Dimensional Problems 

As with two-dimensional solutions, the ft JLLL , term is zero 
due to the cosd term. However, the G mm integral term gives a 
singularity because r is zero in the equation 

je-^/fdxdy. (2.41) 

The singularity can be removed by changing to polar 
coordinates and thus 

(2.42) 


Chertock developed an approach that can easily solve the 
singularity problem in the [G] influence matrix [70]. 
Chertock' s approach finds the radius R 0/ of a circular element 
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that integrates over the same area as a rectangular element. 
Assuming that AR is the area of the m"* element, the singular 
integral becomes 

G =4-(e"* J3nTO -1) • (2.43) 

mm 2 k 

Thus, the singularity in the three-dimensional boundary 
element equation is easily removed. 
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Chapter 3 


SOUND AND STRUCTURE INTERACTION 


Over the past few years, the study of sound- structure 
interaction has received increased attention. The aircraft 
and automobile industries have rigorously studied sound- 
structure interaction effects. Currently the automobile 
industry is pursuing methods to model NVH (Noise, Vibration 
and Harshness) environments. Some components that have been 
modelled in the automotive industry include: carburetor 
manifolds, transmission housings, and tire/hub systems. The 
most popular and most effective coupling method employed for 
the study of their problems is the coupled boundary and finite 
element methods. In industry, most boundary element codes are 
written in-house, while most finite element codes are 
commercial programs. The most common finite element programs 
in use today are NASTRAN, ABAQUS , IDEAS, ANSYS, and ALGOR. 
SYSNOISE was the first program that coupled basic finite and 
boundary element analyses. In 1994, the COMET/BEA programming 
system was introduced into industry. This program couples the 
finite and boundary element methods in an attempt to model NVH 
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environments. Currently, composite laminates and random 
excitation can not be analyzed using the COMET /BEA system. 

This chapter introduces the governing equations for the 
coupled boundary and finite element methods. Coupled BEM-FEM 
results are then compared with existing experimental data for 
a well-known cavity-backed plate problem in Chapter 4. 
Coupled BEM-FEM formulation for plates of composite materials 
is introduced and transmission loss characteristics are 
studied. Random response analysis of the structural and 
acoustic domains is also presented. 

3.1 Acoustic Formulation (BEM) 

In order to couple the boundary element method with the 
finite element method, the governing equations for the 
acoustic domain are slightly modified. The governing 
equations, representing the acoustic domain, for two and 
three-dimensions are the same. Only the influence matrices, 
in Equation (2.31), change as a result of the dimensionality 
of the problem. In this study, only three-dimensional systems 
are coupled. The acoustic domain is governed by the wave 
equation as shown in Equation (2.4). By assuming the time 
dependance as e' ul , the wave equation reduces to the Helmholtz 
equation as seen in Equation (2.7). On applying the boundary 
element method over the domain, the governing equation in 
matrix form becomes 
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(3.1) 


[H]{p> = [G]{|£}. 

The elements of the influence matrices, [H] and [G], are 
given in Equation (2.28). The [H] and [G] matrices represent 
the transfer function of the acoustics and they are functions 
of the input frequency, a>, and the distance vector, r. At the 
acoustic boundary, the normal derivative of pressure is 
related to the acoustic particle acceleration, velocity and 
displacement by 



d 2 u ( 

IP 


du t 

~Zt 


= -ip a w^ T z=p a w i u ( 


(3.2) 


where u, is the acoustic particle displacement at the 
boundary. This boundary condition can be easily derived using 
the continuity equation. On substituting Equation (3.2) into 
Equation (3.1), we get 

[H] {p} = p 0 u 2 [G]{u a }. (3.3) 

Equation (3.3) represents the boundary element method 
governing equation for the acoustic domain. Equation (3.3) is 
applied only on the acoustic boundary of the duct system and 
thus, Equation (3.3) can be modified to couple with the finite 
element method. 


3.2 Structure Formulation (FEM) 

The finite element method has become a very powerful tool 
in analyzing static and dynamic response of structures. The 
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finite element method is capable of handling complex 
structural analysis. The objective here is to model the plate 
or structure and obtain pertinent linear theory information 
needed to couple with the boundary element method. 

3.2.1 Finite Element Specifications 

Many rectangular and triangular type finite elements are 
currently being used in commercial and in-house codes. Any 
type of finite element can be applied to the following 
formulation. The element selected for this study is a four- 
node rectangular element and is shown in Figure 3.1 [67], The 
characteristics of the element are given in Appendix B. The 
element displacement functions for the transverse and in— plane 
directions are 




(3.4) 


and 


u = [H „]{/?} 

v=[ff v ]{0}. (3.5) 

The shape functions, [H.], [HJ and [HJ , and the generalized 
coordinates, a and /?, for the rectangular element are also 
shown in Appendix B. The generalized coordinates are related 
to the element bending and membrane nodal displacements 
through 
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Figure 3.1 Finite Element Rectangular Plate Element. 



(3.6) 


{<*} = [!•,] {urj 

m = lT m ){w m y, 

where the transformation matrices, [T m ] and [T b ] , are functions 
of the element lengths, al and bl. The transformation 
matrices are given in Appendix B along with the definition of 
the degrees of freedom in bending and membrane for the plate 
element. 


3.2.2 Strain-Displacement and Stress-Strain Relations 

The linear relationship between the strain and 
displacement is given as 

{ ® }={£„} + z 0 {k} , (3.7) 

where {e m } is the membrane strain, {*} is the curvature, and 
z 0 is the distance from the mid— surface. The curvature and 
the membrane strain are given by 


and 



- W 'xx 




{*} =< 

~ W ’y> 

► _ 

’rr 

{*} = [C b ]{ a), 


r 

1 

ro 

* 

$ 

\ 


* 

CM 

1 


(3.8) 


(0 = 




U 'y + V ’x 


K 

*v, 

H.. +H 






(3.9) 


where the subscripts b and m denote bending and membrane, 
respectively. The matrices, [C b ] and [C m ], are functions of 
x and y and are given in Appendix B. 
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The linear stress-strain relations for a composite layer 


are given by 




Qn 

£|2 

0 


£ 1 

*2 


Q 2. 

C?22 

0 

* 

C 2 * 

T 12 


0 

0 

C>66 


7,2 


(3.10) 


where the subscript 1 represents the direction of the fibers 
and subscript 2 represents the direction normal to the axis of 
the fibers. The Qy, (i ; j=l,2,6) / components are calculated 
with E u , E 22 , G 12 , v l2 and v 2l . If the lamination angle, <p , is 
changed, the stress-strain relations for the kth layer becomes 


a 


Q\\ Q\2 Q 16 


e 

X 


_ __ 


X 

a > 

► _ 

O 2 I ^22 £*26 

* 


> . 

k 

061 (?62 (?66 

k 



where [Q] is the transformed stiffness matrix. Derivation of 
[Q] and [Q] is given in detail in Appendix B. The stress 
components of the kth layer are given by 


{o} k =[Q) k ({e m } + z 0 {K}) , 


(3.12) 


which are substituted into the stress resultant expression, 


h 

7 

(W,{M fc })= f {a} k (l,z 0 )dz 0 . 

h 

'1 


(3.13) 


This produces the constitutive relations for the forces and 
moments of a composite laminate as 
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(3.14) 




where the [A] , [B] , and [D] laminate stiffness matrices are 

developed from 


h 

1 

( [A] ,[B] ,[£>]) = [ [Q] k (l,z o ,zl)dz 0 . (3.15) 

•'a 

"1 


3.2.3 Derivation of Equations of Motion 

The equations of motion of the plate structure are 
derived using the principle of virtual work. The principle of 
virtual work states that the net work done by a system in an 
equilibrium state under a virtual displacement tends to zero. 
The virtual work definition 
is 


6W=6W inl -SW oa =0 . (3.16) 

The virtual work of the internal forces is 

^=f((<5e m } 7 '{W} + {^> r {M il })dA, (3.17) 

the forces {N} , and the moments {M,,} , are given in 
Equation (3.14). 
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The total external virtual work is 
SW at =^-ph(uSu+vSv+w b Sw b )dA+ ^p(x,y,t) Sw b dA, (3.18) 

A A 


where S u, 6v, and 6w b are the infinitesimal virtual 
displacements. The pressure loading, p(x,y,t) , in Equation 
(3.18) represents the combination of the external pressure and 
the acoustic pressure. The virtual strain and virtual 
curvature in Equation (3.17) are given as 


{6e m }=6{[C m ]{(3}} 
= [C m ]{<S/3}. 


and 


{5k} =S{ [C b ] {a} } 
= [CJ{5a}. 


(3.20) 


Substituting Equations (3.14), (3.19), and (3.20) into 

Equation (3.17), we get 


^=|({^} r [C m ] r ([A]{£ m } + [B]{/c}) 


(3.21) 


+ {*a} r [CJ r ([B]{e m } + [D]{/c})) dA. 

Substituting Equations (3.8) and (3.9) into Equation (3.21), 
we get 
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(3.22) 


Sw ia ,= \({60} T [C m ) T [A][C m U0} 

J A 

+ {60} T [C m ] T [BUC b ]{ a} 

+ {6cc} T [C b ] T [B][C m ]{0} 

+ {Soc} T [C b ] T [D][C b ]{a})dA. 

Equation (3.22) can be written in terms of the membrane and 
bending deflections, and the linear stiffness terms as 


+i*".} T lk,+Hv k } + {6w h } T lk tm }{w m }. 


(3.23) 


The linear membrane and bending stiffness matrices in Equation 
(3.23) are given as 


t*J=[r m ] r [[C m ] r [A] [C m ]dA[T m ], 

u 


and 


[^J=[rj 7 'f[C i ] 7 '[D] [ C b ]dA[T b ] . (3.25) 

The linear coupled stiffness matrices are given by 

(3.26) 

j A 

The in-plane and transverse contributions to the external 
virtual work is given by Equation (3.18). Replacing u, v, and 
Wb with their approximations, and the external virtual work 
can be expressed as 
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M*=-\(<S” m > T lT m VlH u ] r ph[H u ][T m }{w m } 

+ {Sw m } T [T m ] T [H v ] T ph[H v ) [T m ] {tf m } 

+ {Sw b } T [T b ] T [H w ] T ph[H w ][T b ]{w b })dA 
+ {Sw b } T {p b (t) } 

= -{5v m } 7 '( [m^] +[^3 ) {w b } + {6w b } T {p b (t) } . 

(3.27) 

Substituting the expressions for the external and internal 
work. Equations (3.23) and (3.27), into the principle of 
virtual work equilibrium relation, Equation (3.16), we obtain 
the elemental equations of motion 

[m] {w} + [k e/ ] {w} ={p A (t) } . (3.28) 

The linear element mass and stiffness matrices from Equation 
(3.28) are represented as 


and 

[**] = 

The load vector, p b (t) , represents the bending load as 

(P*(t) } = m r jp(x,y,t) {HJdA. (3.31) 


k b k bm 
*mb k m 


(3.30) 


[«] = 


m b 0 
0 m_ 


(3.29) 
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3.2.4 Global Equations and Solution Procedure 

Applying the elemental equation of motion for each 
element and summing each element's contribution, we obtain the 
global equation of motion as 

[M]{tf} + [Jn{tf}={F r }. (3.32) 

Assuming that the time dependance of the plate response is e""*, 
the free vibration problem can be reduced to 

[*]{*>.=« 1 I 2 [*H*} 1I . (3.33) 

Equation (3.33) can now be solved for the natural frequencies 
oj n and the normal mode shapes $ n for the structure. These 
frequencies and mode shapes are used to define the structural 
response function. The response function is then used to 
couple with the boundary element method. 

The damping effects are now added to the structural 
system as 

[Af] + {IV} + [IT] {W}={F t } . (3.34) 

The damping matrix [C] can be experimentally determined or 
assumed as proportional damping. 

A modal transformation is now applied with 

< 3 - 35 > 

where {d> are the modal displacements and [$] is the modal 
matrix. This modal reduction allows the system to be reduced 
to a fixed number of modes. All analyses in this study was 


42 



done using ten inodes for the modal reduction. Ten modes were 
selected to make sure all dominant modes were analyzed. The 
structural equation of motion, Equation (3.34), can now be 
expressed in modal coordinates as 


w of w * ^ <*>" W 


(3.36) 


where M„ and are the modal mass and the damping ratio, 
respectively. From Equation (3.36), a response function is 
easily derived as 




M n (w;;-w 2 +2i $>„(*)) 


(3.37) 

Using the structural response function, the structural 
governing equation becomes 

[HJ {d} = [$] r {F r }, (3.38) 

where the subscript s represents the structural response. The 
subscript used in Equation (3.37) is to distinguish the 
structural response from the influence matrix in Equation 
(2.31). The modal load vector can be separated into external 
forces and the forces applied by the acoustic pressure. 
Therefore, the governing equation becomes 

[H s ]{d} = [i] T {F p } + l*] T {F e } , (3.39) 

where {F p } is the acoustic force and {F e } is the external or 
input force. The acoustic force vector [^] T {F p > , is 
transformed from the boundary element system to the finite 
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element system as shown in Figure 3.2. A transformation 
matrix [LT], is created from the [$] matrix, and is used to 
transfer the acoustic loads from the center of the boundary 
element to the nodes of the finite element or vice versa. The 
[LT] matrix arises due to the difference in the definition of 
the element types. The boundary element is represented by a 
node in the center of the element, whereas the finite element 
is represented by its corner nodes. The [LT] matrix performs 
two tasks simultaneously; it performs the modal reduction and 
the element transformation. Rewriting Equation (3.39) with 
the load transformation matrix, we get 

[H s ] {d} = {LT] T {F p } + [*] T {F e } . (3.40) 

Equation (3.40) is the final form of the governing equation 
for the linear finite element method for the structure. 

3.3 Coupling System Equations and Solution Procedures 
For Harmonic Response 

The structure and acoustic domains have been modelled 
using finite and boundary elements. Equation (3.3) is the 
governing boundary element equation for the acoustic domain 
and Equation (3.40) represents the finite element formulation 
of the structure. A coupled plate/cavity problem, is shown in 
Figure 3.3 to assist in the physical interpretation. The 
acoustic force F p and the external force F c are also shown in 
Figure 3.3. The acoustic domain and the structural domain are 
spatially separated. Coupling the two initially independent 
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systems takes place in terms of compatibility conditions 
across the interface of the plate and the acoustic domain. 
When the compatibility of displacements and forces applied at 
the interface are enforced, the two systems become coupled. 

Applying the displacement compatibility over the interior 
plate surface requires that the acoustic displacement equals 
the out-of-plane plate displacement as 

{W b } = {u a } = [LT]{d}, (3.41) 

which has been related to the modal displacement {d}. The 
definition relation of the out-of-plane displacement is seen 
as , 

<3 ' 42) 

The force compatibility equation represents the force from the 
acoustic pressure as an additional force {F p }, as seen in 
Equation (3.40). Assuming that the acoustic pressure is 
uniform over each plate element, the acoustic force vector, 
from Equation (3.40), may be represented as 

[LT] T {F p } = [LT] T [S A ){ Pacl } , (3.43) 

where {p, cl } is the global acoustic pressure vector represented 
at the center of the boundary element. The diagonal matrix 
[S A ], is composed of the elemental areas of each interface 
element . 


47 



The external force vector from Equation (3.40) may also 
be expressed in terms of pressure as 

[«] r {^} = -[$] r {P t (t)}=-p w [$] r {P M (t)}=-p w {f mod }, 

(3.44) 

where {p^} is the constant input pressure vector and {f^} 
represents the modal force vector. The vector (P b (t)}, is the 
global finite element force vector, which is composed of the 
contribution of the element forces p b (t) , from Equation 
(3.27). The value p^ has been pulled from {P b (t)> in Equation 

(3.44) , since uniform pressure loads are assumed. In all 

transmission loss calculations, a unit input pressure vector 
is assumed. The acoustic and structure system can now be 
coupled, using Equations (3.3), (3.40), (3.41), (3.43) and 

(3.44) . In matrix form, the coupled system is given as 


LH] 

-p a o) 2 [G] [LT] 


r 

Pact 


0 

-[L!T] r [SJ 

[*,] 


d 


od> 


Equation (3.45) can now be solved for all acoustic pressures 
and modal displacements. The global pressure vector {p act }, 
encompasses all pressures from the entire acoustic boundary. 
For instance, if the acoustic boundary is discretized into 150 
boundary elements and the plate into 25 finite elements, the 
pressure vector {p, cl }, contains all 150 acoustic pressures. 
Only the first 25 pressures are coupled with the plate system. 
Note that this system is only a function of frequency. After 
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solving for the acoustic boundary pressures and the modal 
displacements, a variety of other data becomes readily 
available. From Equation (3.35), the element plate velocities 
can easily be calculated. Using the boundary pressures {p lct } 
as boundary conditions, any interior acoustic pressures or 
velocities can be calculated from Equation (3.3). 

The importance of coupling the acoustic and structural 
domain will be shown in Chapter 4. In order to accurately 
model aerospace and automotive vehicle environments, coupled 
systems must be considered. This study is a step towards 
improving the modelling of such complicated environments. 

3.4 Coupling System Equations and Solution Procedures 
For Random Response 

The majority of real world systems involve random loading 
and random response. As stated in the introduction, this 
study is an attempt to accurately model structures in complex 
environments for aerospace and automobile industries. In 
order to satisfy this criteria, the previous harmonic analysis 
is now modified to include random loading and response. The 
goal of this section is to analyze the previously selected 
problem and model it with linear random analysis. The random 
response of the coupled plate/cavity system will be compared 
with the random response of the plate structure only. 
Therefore, a brief introduction of the random response of 
composite panels is now given. 
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3.4.1 Random Response of Plate Structure 

The governing finite element equation for the linear 
plate structure is given in Equation (3.34). A standard 
eigenvalue/eigenvector solver is applied with the global mass 
and stiffness matrices and the natural frequencies and mode 
shapes can be calculated. The modal equation has been given 
in Equation (3.36). The spectral density of the plate 
response is given by 

(S*(«) ),= !*,(«> \ 2 n S w (w) , (3.46) 


where 


W«>>.- r-ri 

M n (u„- <a l +2if m <u H <a) 


(3.47) 


The s m (u) value is the input spectral density of the external 
excitation P b , from Equation (3.44). Equation (3.46) holds 
only for linear random analysis. From the definition of 
spectral density, the mean square response of the modal 
displacement: can be written as 


or 


E(dl).S p {u) (. f moi ) 2 n 


jj/Uco) \l d<a, 


(3.48) 


E{d\) 


**, (<■»(*■ Cd) 

4 LMWn 


(3.49) 
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The input spectral density is assumed as white noise and 

S p ( w) =8.4216 X 1 ->•“», (3.50) 

where SSL is the input sound pressure level in dB. The mean 
square of the modal displacement can be related to the mean 
square of the displacement by 

E[{^} {W} r ] =[*]£[ {d} {d} r ] [*] r , (3.51) 

where [$] is the modal matrix [67]. Since the terms E[d„dj] 
(n?*j)are at least two orders of magnitude smaller than Efd,, 2 ] 
(n=j ) , all cross-correlation terms, Etc^dj] (n*j),are set to 
zero. Therefore, Equation (3.51) reduces to 

Hof modes 

E(wl) - £ (3 - 52) 

j - 1 

From Equation (3.52) , the maximum root mean square (RMS) 
deflection can be calculated and the results can be compared 
with those from the coupled analysis results. A detailed 
derivation of the random response of plate structures can be 
found in [67] and [68]. 

3.4.2 Random Response of Coupled Structural-Acoustic 
System 

The random response of the acoustic domain is formulated 
from the wave equation [68]. A classical continuum approach 
is used to relate the acoustic pressure and the modal 
displacement through a transfer function. The fundamental 
pressure response of the acoustic domain can be written as 
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(3.53) 


P(x,y,z,t) =H p ^ / {x ,y ,z ,u)u a (t) , 

where Hp^ is the acoustic transfer function. In other words, 
the acoustic pressure at any location in the duct can be 
related to the acoustic displacement through the acoustic 
transfer function. In order to find the acoustic transfer 
function, u, is set to le‘"\ Substituting Equation (3.53) into 
the three-dimensional wave equation, we get 

^ 2 Hpact + ^Hpaa = 0, (3.54) 

where k is the acoustic wavenumber and u, represents a unit 
harmonic response. The boundary conditions for Equation 
(3.54) for the coupled system in Figure (3.3) are 


tepaa 

dn 


= 0 


on a rigid wall 


8H 


pact __ 


dn 


_ 2 „ for a displacement 

Pa “ at the interface . 


(3.55) 


To find the transfer function Hp,,.,, the displacement vector is 
set to unity. Equation (3.54) can be applied to the acoustic 
domain utilizing the boundary element method, as 


dH 


(3.56) 


An expression for the acoustic transfer function can be 
obtained by 
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(3.57) 


{ V«} = W '* ^ C LT 3 < "^“ 2 > =< C c> » 

where the [G] matrix has been modified to account for only the 
vibrating wall of the duct and the element transformation has 
taken place through Equation (3.41). 

By applying the boundary element method, a set of 
equations is created where {Hp, ct } is now only a function of 
frequency and not space (x, y and z) . The boundary conditions 
for the acoustic domain. Equation (3.55), have been applied to 
Equation (3.56), which yields a solution for the acoustic 
transfer function, {H„, ct } . The vector {Hp« t } represents the 
acoustic transfer function for each boundary element due to a 
unit displacement at the structure. For linear random 
analysis, the spectral density of the pressure response can be 
related to the input spectral density by 

{^(w)} = {H pacl (u)}{H pacl {-u)} T {S d {u)} {3 58) 

= {C c (u)} {C c (-a)} T {S d (a)} . 

Equation (3.58) is used to model the random response of the 
acoustic domain and it relates the spectral density of the 
acoustic pressure to the spectral density of the input 
acoustic displacement. 

The structural system can be modelled through Equation 
(3.46), taking into account the acoustic spectral density and 
the input force or external force spectral density. Assuming 
that the input force spectral density is a constant white 
noise input and that the spectral density of the acoustic 
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pressure is also constant over each element, Equation (3.46) 
reduces to 


(S d (u)) n =\H s (u) \ 2 n S p (u>) (f^l 

, « , (3.59) 

-!*,(«> I- E ( S poc(«))y/ 

where NPEL is the number of plate elements. In order to 
transfer Equation (3.59) to matrix form, we write 

M of modes NPEL 

[LOAD] = £ £ ( LTj n ) . (3.60) 

n«l jm 1 

The coupled system can now be modelled using Equations (3.58) , 
(3.59) and (3.60) as 


I 

[LOAD] |H,(oj) 


-{C c (u) }{C f (-o>)} r 

I 


{^(co) | - 

^(«){^od}|^(«) I'}' 


(3.61) 

where [I] is the identity matrix. Equation (3.61) couples the 
spectral densities of the plate\cavity system in the frequency 
domain. The unknowns in Equation (3.61) are the acoustic 
pressure spectral densities, {S pic J, for each acoustic boundary 
element and the structure modal displacement spectral 
densities, {S d }. 
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For simplicity, Equation (3.61) can be written as 

[ CPLG(u )] {SD(w)> = {SD ;JV (o>)}, (3.62) 

where SD(co) represents the unknown spectral densities from 
Equation (3.61). The matrix [CPLG(u)] represents the total 
coupling matrix in Equation (3.61). Multiplying both sides of 
Equation (3.62) by the inverse of [CPLG] and integrating over 
the derived frequency range, (-oj c <co<oj c ) , we get 

u r 

(SD(w) }du>= | [CPLG(w) y l {SD m (u) }du > . 

-CJ f 

(3.63) 

The integral in Equation (3.63) can be solved using a simple 
trapezoidal approximation. After applying the integration in 
Equation (3.63), the mean square response of the modal 
displacement and the acoustic pressure can be readily 
calculated. 
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Chapter 4 


NUMERICAL RESULTS 

In this chapter we present the results of the previous 
formulation and also provide some additional discussion. The 
results section has been subdivided into three sub-sections: 
acoustic results, coupled BEM-FEM results, and analysis of 
panels in the Thermal Acoustic Fatigue Apparatus (TAFA) at 
NASA Langley Research Center. The purpose of this chapter is 
to show that the boundary element method is a powerful 
numerical technique that can be used to solve a wide range of 
acoustic problems. The coupling of the boundary element and 
finite element methods, leads to a powerful tool that can be 
readily used to analyze structural-acoustic problems. 

4.1 Acoustic Duct Results 

The boundary element method has become a preferred method 
when modelling acoustic domains. BEM has advantages such as: 
reduction of discretized dimensionality of the domain, 
implementation of infinite domains, and shorter computation 
time. Here we will analyze two and three-dimensional 
interior acoustic domains. 
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4.1.1 Two-Dimensional Interior Domain Analysis 

Due to the symmetry of structures, a two-dimensional 
analysis can be performed to save time and calculations. Very 
little work has been done in the study of 2-D acoustic fields 
using the boundary element method due to difficulties 
associated with singular Bessel functions. Therefore, some 
acoustic field problems with exact solutions are selected to 
see how well the boundary element method can approximate these 
solutions. Some very useful information is learned about the 
number of elements required to solve problems and the limits 
on the wavenumber for accurate results. 

The first problem is a simple rectangular duct with an 
open outlet. The geometry is shown in Figure 4.1, and the 
boundary conditions are given as follows: 

|£=0 at y=0 and lm 

, ^ „ (4.1) 

p=l at x=Om 

p = 0 at x=6 m. 

The duct length is set at 6 meters and the wavenumber k, is 
1. The exact solution is 


p(y) = sin k (L-x) 
sin kL 


/ 


(4.2) 

where the pressure p(x) only depends on the x coordinate since 
a plane wave solution is assumed. The pressure along the 
boundary is calculated versus the x-distance. The comparison 
between constant and linear boundary elements with the exact 
solution is shown in Figure 4.2. Eighty boundary elements are 
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used for both the linear and constant element cases. The 
element distribution is as follows: 32 boundary elements on 

y=0.0, 32 boundary elements on y=1.0 m. , 8 boundary elements 
on x=0.0 and 8 boundary elements on x=6.0 m. All elements are 
of equal length on each individual surface. As expected, the 
linear elements give a better approximation to the exact 
solution than the constant elements. 

Also for the open outlet case, the pressure is plotted at 
(x,y)=(l. 5, 0) versus the wavenumber. The wavenumber is varied 
from 0 to 4 as shown in Figure 4.3. The natural frequencies 
of the cavity appear at the pressure discontinuities. Figure 
4.3 shows the pressure variation for linear elements only. 
The exact cavity frequencies can also be easily calculated. 
The acoustic natural frequencies for constant and linear 
elements are compared with the exact frequencies in Table 4.1. 

the constant elements and the linear elements give very 
good approximations to the acoustic natural frequencies. 
Knowing that the boundary element method can approximate the 
accurately at the acoustic frequencies becomes very 
important when coupling BEM and FEM. 

The convergence of the pressure values versus the number 
of elements is also very important. The same type of problem 
is considered except that the length of the duct in the y 
direction is set at 6 meters which is the same as the x 
direction length. The wave number is set at .001 to ensure 
the plane wave assumption and the pressure was observed at 
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(x,y)=(3,5) for linear elements. The convergence data is 
shown in Table 4.2. As the number of elements increases, the 
error sharply decreases as expected. 

The second example problem used to validate the boundary 
element method is a rectangular duct with a closed outlet as 
shown in Figure 4.4. The closed outlet is assumed to be a 
rigid wall at the exit. The boundary conditions for the duct 
in Figure 4.4 are, 


^£=0 at y=0 , lm 
an 

^£=0 at x=6 m 
an 


(4.3) 


p=l at x=0m . 

As before, the wavenumber is set at 1 and the length of the 
duct in the x-direction is 6 meters. The exact solution is 
easily found as, 


CCS k(L-x) . (4.4) 
cos kL 

The pressure, p(x), is only a function of the x-direction due 
to the assumption of plane wave propagation. As before, 80 
elements are used for the constant and linear element cases. 
The boundary element comparison with the exact solution is 
shown in Figure 4.5. Both constant and linear element cases 
give very good approximations. 

In most realistic duct problems, the velocity or the 
pressure is not known at the duct outlet. Only the 
relationship between the normal pressure and the normal 


61 



Table 4 . 2 


Pressure Convergence with Number of Elements. 


Number of Linear Elements 

Pressure 

4 

0.47561 

8 

0.49979 

12 

0.49987 

80 

0.50000 

Exact 

0.50000 




is known. This is known as an impedance boundary 
condition, where z=p/v, is some arbitrary complex constant. 
A realistic condition at the exit requires that no reflected 
acoustic waves are encountered at the outlet. Since a plane 
wsve assumption is made, the no— ref lection exit condition 
requires that the impedance be set to p,c, where p t is the 
®®dium density and c is the medium sound speed. Knowinq the 
relationship between the pressure and the velocity at the 
ex it , 3 boundary condition in terms of p and q can be derived. 

The next example problem includes an impedance boundary 
condition at the outlet. The example problem in Figure 4.6 
includes the following boundary conditions: 


3p_ 


IE 


=0 at y= 0, lm 


z = — =P„c at x=6m 

V ' “ 


p=l at x=0m. 


(4.5) 


The duct is discretized into 80 elements for both constant and 
linear element cases and the wavenumber is set at 0.05. In 
Figure 4.7, the pressure results for the linear and constant 
element cases are plotted along with the exact solution versus 
the x distance of the duct. Both element cases give extremely 
accurate results. The constant element case gives 
continuously higher values than the linear case, however the 
maximum error for the constant element case is only 0.3 
percent. The results of this case are very good because the 
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Figure 4.7 Results for a Rectangular Duct with an 
Impedance Outlet Condition. 
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acoustic wavelength is large compared with the duct element 
dimensions. 

The boundary element method has shown that it can handle 
any type of boundary condition of pressure and velocity on the 
surface of the duct system. The next objective is to see how 
BEM handles sound propagation in rectangular ducts with 
complex geometries. A muffler— type problem with an exact 
solution is modelled and compared with the approximations of 
the boundary element method. 

The variable geometry of the muffler-shaped duct is shown 
in Figure 4.8. The side walls are considered rigid and there 
is a unit pressure distribution at the inlet of the duct. In 
order to assume no reflection at the outlet, an impedance 
boundary condition is used, where z=p/v 1 =p,c. The 
transmission loss coefficient is calculated between the input 
pressure and the outlet pressure. The cross-sectional areas 
S, and S 2 , are set at 1.0 m and 2.0 m respectively. Thus the 
ratio of the two areas, S 2 /S,, is set at 2.0 and the 
rectangular duct is discretized into 80 elements. The inlet 
section and the outlet section use 20 boundary elements each 
and the mid-section of the duct uses 40 boundary elements for 
the discretization. The length of section 2 as shown in 
Figure 4.8 is set at 1^=1. 0 m. 

The exact solution for the transmission coefficient can 
easily be derived in terms of complex amplitudes of incident, 
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reflected, and transmitted pressures [64]. At low frequencies 
corresponding to kL^l, leads to a transmission coefficient of 


“. = ■ 


1 + 


S 7 kL 7 

2S, 


(4.6) 


Equation (4.6) shows that the transmission coefficient is near 
unity for low frequencies and tends to zero as the frequency 
is increased. The equations derived for pressure transmission 
are valid only when the wavelength is large compared with the 
smallest dimension of the duct section. 

As the wavenumber is increased, the transmission 
coefficient is evaluated for both constant and linear elements 
as shown in Table 4.3. The transmission coefficient decreases 
as the frequency increases as expected. Both constant and 
linear elements give very good approximations with average 
errors of 1.36% and 1.08%, respectively. The error values are 
calculated up to k=.4 because the exact solution may not be 
valid due to the kLj« 1 assumption for larger k values. 

From the pressure comparisons in the above example 
problems, it can be seen that the boundary element method is 
a very good approximating method for two-dimensional domains. 
The next objective is to show the versatility of the boundary 
element method through an example. A multi-section duct was 
created as shown in Figure 4.9. The side walls are considered 
to be rigid and a unit pressure field is applied at the 
entrance of the duct. An impedance boundary condition is 
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Figure 4.8 Muffler-Shaped Duct System. 


Table 

4.3 Pressure 

Values for Muffler- 

-Shaped Duct. 

— 

k 

Exact Sol. 

Linear Elem. 

Constant Elem. 


0.01 

0.9999 

0.9999 

0.9999 


0.1 

0.9901 

0.9959 

0.9965 


0.2 

0.9615 

0.9806 

0.9829 


0.3 

0.9174 

0.9437 

0.9475 


0.4 

0.8621 

0.8622 

0.8679 

- 

0.5 

0.8000 

0.6941 

0.6985 

— 
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employed at the outlet to assume no acoustic reflection and 
the wavenumber is set at 0.1. The three-duct system is 
discretized into 70 boundary elements and the transmitted 
pressure is observed at the center of the outlet duct. 

The transmitted pressure is calculated for both constant 
and linear elements as the angle <p l is varied from 0 to 45 
degrees. The results are plotted in Figure 4.10 for angle 
values, 0,, in increments of 5 degrees. Both constant and 
linear element cases give very close results due to the low 
frequency assumption applied at the beginning. In Table 4.4, 
a frequency-pressure response is tabulated for different 
values of the wavenumber. The angle of the mid-section duct 
is set at 20 degrees and the transmitted pressure is obtained 
at the center of the outlet section. The transmission loss is 
calculated with the relationship between the input pressure 
and the output transmitted pressure, TL=201og (Pb/Po,,,) • A 
negative transmission loss value means the output pressure is 
greater than the input pressure. The distance of the mid- 
section duct in Figure 4.9 is set at 1.0 m. and the ratio of 
cross-sectional areas S 2 /S,, equals 0.25. As the wavenumber 
tends towards unity, higher order modes begin to propagate and 
the transmitted pressure begins to rapidly change. Due to the 
odd geometrical configuration of the duct system, higher 
order modes begin to propagate at lower frequencies than 
typical rectangular duct systems. 
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Table 4.4 Pressure Values at the Center of the Outlet 
Section of a Three-Section Duct with 0,=2O°. 


k 

Constant Elements 

Linear Elements 

0.001 

0 .24735 

0.22416 

0.01 

0.24736 

0.22419 

0.1 

0.25617 

0.23268 

0.2 

0.28206 

0.25893 

0.3 

0.29800 

0.29063 

1.0 

-0. 18069 

-0 . 18105 
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As seen from the comparisons above, the boundary element 
method gives very accurate and reliable solutions for the two- 
dimensional duct analysis. Bessel function singularities have 
been approximated. 

4.1.2 Three-Dimensional Interior Domain Analysis 

Acoustic analysis on most realistic structural components 
are three-dimensional in nature. Since the components of the 
wavenumber analysis are affected by the dimensionality, three- 
dimensional analysis must be performed. For example, the 
majority of commercial mufflers, particularly the ones 
embodying flow-reversing chambers, are three-dimensional. 
Thus, one and two-dimensional frequency-domain acoustic 
theories as well as time-domain finite wave theories are not 
applicable [65]. One has to resort to the numerical solution 
techniques in order to solve such complex systems. 

As mentioned before, three-dimensional acoustic analysis 
is easier to apply than the two-dimensional analysis simply 
because the approximate function deals with exponential 
functions instead of Bessel functions. This section first 
shows that the boundary element method is accurate by 
comparing it to some problems with exact solutions. Next some 
example problems with complex geometry features are analyzed 
using BEM to show the versatility of the method. 
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The example problems are derived from a typical 
rectangular duct as shown in Figure 4.11. The duct shows how 
the boundary is discretized with two-dimensional elements. 

Each example problem assumes that each side of the duct 
is 2 0 cm and the entire duct is discretized using boundary 
elements. The two end surfaces (inlet and outlet) are each 
modelled by a 5x5 mesh or 25 elements, and the duct side walls 
are also modelled by a 5x5 mesh or 25 elements each. Thus, 
the total number of elements is 150 constant boundary 
elements, as shown in Figure 4.11. The frequency is varied in 
each case and the pressure is measured at the center of the 
duct outlet. A unit pressure field is assumed at the duct 
entrance (z=0.0 cm) and all of the exterior walls are assumed 
rigid. For the first example problem, a null pressure field 
is applied at the duct outlet and the transmission loss is 
plotted versus the input frequency in Figure 4.12. 

In each case, the input acoustic waves are considered to 
be plane waves and thus the one-dimensional exact solutions 
are valid. From Figure 4.12, the boundary element method 
gives very accurate results compared to the exact solution. 
However, at the acoustic natural frequency of the duct system, 
the boundary element method results vary from the exact 
solution slightly. The cut-off frequency is selected at 1200 
Hz due to the accuracy of the boundary element method. This 
frequency dependent accuracy will be discussed at the end of 
this section. For three-dimensional domains, only constant 
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FREQUENCY (Hz) 

EXACT SOLUTION + BEM SOLUTION 

Figure 4.12 Results for Three-Dimensional Duct 
with Null Pressure at Outlet. 
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elements are used since the two-dimensional analysis showed 
that the discrepancy between linear and constant elements was 
nominal . 

The second example problem uses the same duct dimensions 
and the same acoustic assumptions. The outlet boundary 
condition assumes a rigid wall at z=20 cm. By assuming a 
rigid wall boundary condition, the normal velocity and q are 
set equal to zero. The results for a rigid wall boundary 
condition are shown in Figure 4.13. Once again, at the 
acoustic natural frequency, the boundary element method gives 
erroneous results. BEM gives very accurate approximations, 
except at the acoustic natural frequencies, because the 
influence matrices become singular. 

The last example problem applies an impedance boundary 
condition at the outlet of the duct. A no-reflection 
condition is used, which assumes the relationship between p 
and v equals p,c for the plane wave assumption. The density 
for air is 1.21 kg/m 3 and the speed of sound in air is given 
as 343 m/s 2 . Obviously, the no-reflection condition assumes 
the total plane acoustic wave is transmitted at the outlet. 
The results for the impedance boundary condition case are 
shown in Figure 4.14. 

Notice that as the frequency increases, the boundary 
element results tend to diverge from the exact solution. This 
is a result of two things. Once again, an increase in the 
number of elements for higher frequencies must be applied to 
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4.13 Results for Three-Dimensional Duct with 
Rigid Wall Outlet. 
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4.14 Results for Three-Dimensional Duct with 
Impedance Boundary Condition. 
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assure accurate results. 


This will be discussed later. 


Second, at higher frequencies, the higher order modes come 
into play in the BEM approximations and thus the comparison 
with the exact solution tends to diverge. Therefore, these 
two phenomena contribute to the divergence of the boundary 
element results. 

From the previous results, the boundary element method 
has shown that it gives accurate results for three-dimensional 
domains. The next objective is to model a three-dimensional 
duct with complex geometry. Munjal has solved an offset 
inlet-outlet duct as shown in Figure 4.15 [64]. The inlet of 
the duct system enters at the bottom right of the main section 
and the outlet exits at the top left. Munjal approximates the 
transmission loss of the system with 40 finite elements. 
Three case studies are considered with the boundary element 
method: 50, 82, and 328 elements, as shown in Figure 4.16. 
The 50 element case discretizes the inlet and outlet sections 
with 5 boundary elements each and the mid-section with 4 0 
boundary elements. The 82 element case discretizes the inlet 
and outlet sections with 9 elements each and the mid-section 
with 64 boundary elements. The 328 element case discretizes 
the inlet and outlet sections with 36 boundary elements each 
and the mid-section with 256 boundary elements. The 50 
element case gives a good general approximation. The 82 and 
328 element cases give extremely close approximations to the 
finite element results. Notice that the 50 element BEM case. 
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Figure 4.15 Three-Dimensional Offset Inlet-Outlet 
Duct. 


Transmission Loss (dB) 



FEM- 40 Elements + BEM- 50 Elements 

* BEM- 82 Elements 0 BEM-328 Elements 

Figure 4.16 Transmission Loss of Offset Inlet-Outlet 
Duct . 
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gives erroneous results at very low frequencies. This 
phenomena has been noted before by Hussain and Peat [22]. 
They proposed a small change in their formulation for low- 
frequency calculations. Careful consideration has to be taken 
into account when modelling complex duct systems. 

With confidence in modelling three-dimensional domains, 
a very complex system can now be solved. In Figure 4.17, the 
complex duct system has four inlets and four outlets. All 
inlets and outlets have the same dimensions and are located at 
the corners of the main duct system. Two boundary element 
cases are considered: transmission loss for 74 and 130 
elements as shown in Figure 4.18. The 74 element case 
discretizes each inlet and outlet section with 5 boundary 
elements each and the mid-section with 24 elements. The 
remaining 10 boundary elements are used to discretize the 
inlet and outlet surfaces. The 130 element case discretizes 
each inlet and outlet section with 9 boundary elements each 
and the mid-section with 48 elements. The remaining 10 
boundary elements are again used to model the inlet and outlet 
surfaces. A unit pressure field is activated at each inlet 
and an impedance boundary condition, z =p,c, is assumed at each 
outlet boundary. At various frequencies, the 74 element case 
gives varying transmission loss results from the 130 element 
case. However, both cases give very similar approximation 
trends . 
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Figure 4.18 Four-Inlet and Outlet Duct Results. 
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From the previous results, the boundary element method 
has shown that it gives very accurate solutions for regular 
and irregular-shaped three-dimensional ducts. Some 
interesting aspects have been noted. Some erroneous results 
can occur with low-element, low-frequency modelling. It has 
been found that four to five elements per acoustic wavelength 
are needed to accurately model acoustic domains. Note that 
this generalization applies only to constant elements. It is 
also shown at certain acoustic natural frequencies, that the 
boundary element method can give erroneous results. Many 
authors have addressed and solved this problem with the 
"CHIEF" and "CONDOR" methods [8,48]. These methods solve the 
singular matrix system and provide very accurate solutions at 
natural frequencies. 

At this point, the boundary element method has shown that 
it gives accurate solutions for two and three-dimensional 
domains. Some very complex geometric systems have been 
modelled. The most attractive advantage of the boundary 
element method is that only the boundary is discretized. This 
cuts down considerably on user modelling time. 

4.2 Coupled Harmonic Response of Structural- 
Acoustic Problems 

Coupling the structure and acoustic domains is essential 
when trying to accurately model harsh environments. 
Previously, isotropic panels have been coupled with acoustic 
domains as mentioned in Chapter 1. The first objective of 
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this section is to validate the coupled BEM-FEM approach by 
comparing with known analytical and experimental results for 
isotropic panels. Two selected problems with different plate 
boundary conditions are used for the validation. 

Composite panels are then coupled with the acoustic 
domain and transmission loss characteristics are studied. No 
exact solutions or experimental data are available for sound- 
structure interaction problems using composite panels. Random 
analysis results are included in the last section. 

4.2.1 Isotropic Structures 

As mentioned in Chapter 1, analytical, experimental and 
exact solutions are available for isotropic panels backed by 
an acoustic cavity as shown in Figure 3.3. Guy and 
Bhattacharya provided modal summation approximations along 
with valuable experimental data for two coupled problems [57]. 
Bokil developed exact solutions using Laplace transforms for 
the same two coupled problems [58]. The data from references 
[57] and [58] are very valuable because they give coupled 
freguencies to compare with the BEM-FEM procedure. The first 
cavity backed structure (see Figure 3.3), is a simply 
supported brass panel with the following plate and duct 


characteristics : 



Panel/Cavity Problem - Simply supported Brass Plate 

Cavity 

X axis length = 20 cm. 

Y axis length = 20 cm. 

Z axis length = 20 cm. 

All walls are considered rigid 
Air medium: p,=1.21 kg/m 3 , c=343 m/s 

Plate 

X axis length = 20 cm. 

Y axis length = 20 cm. 

Thickness = 0.9144 mm. 

Young's Modulus = 10.4E10 Pa. 

Poisson's Ratio = 0.37 

Density = 8500 kg/m 3 . 

Two element-mesh cases were run for this coupled system, 
one case consisted of 25 (5x5 mesh), finite elements and 150 
(5x5 mesh for each of the 6 duct walls), boundary elements, 
and the second case used 64 (8x8 mesh) , finite elements and 
384 (8x8 mesh for each of the 6 duct walls) , boundary 
elements. The pressure was calculated at the back wall of the 
cavity (p,*,,) and the transmission loss was calculated using the 
source pressure as the input pressure (p^) . The coupled BEM- 
FEM results are compared in Figure 4.19 with approximate 
analytical, experimental and exact solutions from references 
[57] and [58]. The results in Figure 4.19 indicate that the 
coupled boundary element/f inite element procedure reproduces 
results that are of comparable accuracy with existing 
analytical, experimental and exact solutions. 

The second case considered analyzes the same plate/cavity 
system as shown in Figure 3.3. However, this time the panel 
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Figure 4.19 BEM-FEM Comparison for Simply Supported 
Sound-Structure Interaction Problem. 
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is made of aluminum and clamped instead of simply supported. 
The cavity and plate dimensions and properties are given 
below: 


Panel/Cavity Problem - Clamped Aluminum Plate 

Cavity 

X axis length = 30.48 cm. 

Y axis length = 15.24 cm. 

Z axis length = 15.24 cm. 

All walls are considered rigid 
Air medium: p a =1.21 kg/m 3 , c=343 m/s 

Plate 

X axis length = 30.48 cm. 

Y axis length = 15.24 cm. 

Thickness = 0.16256 cm. 

Young's Modulus = 7.0E10 Pa. 

Poisson's Ratio = 0.33 

Density = 2400 kg/m 3 . 

The transmission loss is once again plotted in Figure 
4.20 with the analytical, experimental and exact results from 
references [57] and [58]. The same two element-mesh cases 
were used for the boundary and finite elements. Again, both 
element-mesh cases give very good results. 

When studying the coupled system, the coupled natural 
frequencies are of vast importance. In Table 4.5, a 
comparison of the first few natural frequencies between the 
coupled BEM-FEM procedure and references [57] and [58] are 
given. The uncoupled natural frequencies are given, where the 
numbers in bold represent the uncoupled plate frequencies and 
the numbers in italics represent the acoustic uncoupled 
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Table 4.5 Comparison of Natural Frequencies (Hz) 
for Coupled Boundary Element and 
Finite Element Method. 


Modal 

Exact 

Exper . 

BEM-150 

BEM-384 

Uncoupled 


[57] 

[58] 

[57] 

FEM-25 

FEM-64 

Plate/ Acoustic 

— 


Simply Supported Brass 

Plate 



91.5 

87.0 

91.0 

85.0 

86.0 

78.1 

- 

390.0 

390.4 

397.0 

384.0 

385.0 

390.3 


695.0 

702.5 

730.0 

697.0 

705.0 

702.5 


860.0 

860.0 

864.0 

860.0 

861.0 

860.0 

— 



Clamped 

Aluminum Plate 



238 . 1 

234.4 

235.0 

234.0 

238.0 

227.0 


564 . 0 

588.5 

565.0 

578.0 

562.0 

590.1 
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frequencies. The uncoupled natural frequencies are given in 
order to show how the coupling affects the total system 
response. Note that the first natural frequency of the 
coupled system with the simply supported plate is greater than 
the first uncoupled plate natural frequency. This is because 
the acoustic chamber "spring-loads" the plate and thus 
increases the coupled natural frequency. Also note that the 
coupled system is affected by the first acoustic natural 
frequency at 860 Hz. The same phenomena occurs with the first 
coupled natural frequency of the clamped aluminum panel. 
However, note that only the first natural frequency is 
affected by the coupled system. 

From the two systems above, the coupled BEM-FEM procedure 
gives very reliable results. Both element cases seem to 
favorably approximate all of the analytical, experimental and 
exact solutions. 

Another interesting aspect of the structural-acoustic 
system is the effect of the cavity depth on the plate mode 
shape. Figure 4.21 shows the first mode shape of the square 
brass plate driven at the coupled natural frequency of 86 Hz. 
Increasing and decreasing the length of the acoustic chamber 
greatly affects the mode shape of the panel. If the acoustic 
chamber is reduced from 20 cm. to .2 mm., the plate appears to 
vibrate like the third mode shape as shown in Figure 4.22. 

The outlet boundary condition can provide another 
interesting study. In all of the above problems, the outlet 
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Figure 4.21 First Mode Shape of Brass Plate 
Vibrating at 86 Hz with Duct 
Length of 20 cm. 
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Figure 4.22 Third Mode Shape of Brass Plate 
Vibrating at 86 Hz with Duct 
Length of .20 mm. 
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boundary condition is assumed to be a rigid wall. Here, we 
will assume that the outlet boundary condition has a constant 
impedance value. Thus the relationship between pressure and 
velocity is varied from p,c to 1, ( . 0415 < z < 1 . 0) . As z 

increases towards 1, the rigid wall approximation is applied. 
When z equals p t c, there is no impedance mismatch at the 
outlet, and the acoustic wave travels out of the duct system. 
From Figure 4.23, the coupled natural frequency is equal to 
the uncoupled natural frequency of the panel (78 Hz) in the 
low impedance range. As the impedance value increases, the 
coupled natural frequency increases up to the maximum coupled 
natural frequency of 85 Hz (the same as the rigid wall 
condition) . 

4.2.2 Composite Structures 

As the next step towards exploring the area of sound- 
structure interaction, this study offers numerical results on 
an acoustic cavity-backed composite panel by harnessing the 
strengths of the coupled BEM-FEM technique. The composite 
panel used in this study has the following properties: 

Graphite /Epoxy 

E n = 15.57E10 Pa 
E 22 = 0.807E10 Pa 
v \2 ~ 0.22 
G 12 = 0.455E10 Pa 
p = 1550 kg/m 3 . 


Natural Frequency (Hz) 



Figure 4.23 Coupled Natural Frequency Results 

With Variable Impedance at the Outlet. 
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The first transmission loss analysis is performed on the 
coupled system shown in Figure 3.3 and the same dimensions 
(20x20x0.09144 cm) are used from the Brass Panel problem in 
the previous section. The composite panel is simply supported 
and four panels with a different number of symmetrical cross- 
ply layers are considered. The number of layers is varied 
from 1, 3, 5 and 7. The one-layer laminate is a 0° degree 
layer, the three-layer laminate has a (0°/90 o /0°) layup, the 
five-layer laminate has a (0°/90 o /0 o /90 o /0°) layup and the 
seven-layer laminate has a (0 o /90 o /0 o /90 o /0°/90 o /0 0 ) layup. 
Figure 4.24 shows the transmission loss comparison for the 1 
and 7-layer composite panels and the brass panel. The 
transmission loss data for the 3 and 5-layer composite panels 
lay between the curves of the 1 and 7-layer laminates and thus 
are not shown. The first coupled natural frequency of the 
composite panels is much greater than the coupled natural 
frequency of the brass panel. Note that the second coupled 
natural frequency for the one-layer composite is less and for 
the seven-layer composite is greater than the second coupled 
natural frequency of the brass panel. On an average, the 
brass panel give a greater transmission loss than the 
composite panel. However, remember that the brass panel is 
about 5% times heavier than the composite panels. 

The composite panels are also utilized in the second 
problem of the previous section and transmission loss data are 
compared with the aluminum clamped panel in Figure 4.25. The 
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Figure 4.24 Composite Panel Analysis of Simply 

Supported Structure-Acoustic Interaction 
Problem. 
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Transmission Loss (dB) 
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Composite Panel 

Figure 4.25 Composite Panel Analysis of Clamped 

Structure-Acoustic Interaction Problem. 
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transmission loss of one-layer and seven-layer composite 
panels are again shown along with the aluminum panel curve. 
As before, the three-layer and the five-layer composite panel 
transmission loss curves are bounded by the 1 and 7-layer 
curves and thus are not shown. Note, for this clamped case, 
the first coupled natural frequency of the composite laminates 
is less than the coupled natural frequency of the aluminum 
panel. The second coupled natural frequency of the composite 
laminates is greater than that of the aluminum panel. This 
effect occurs because the first coupled frequency is 
considered "mass-loaded" and the second frequency is 
considered "stiffness-loaded". From Figure 4.25, the 
composite panels give approximately the same transmission loss 
characteristics in an average sense. Remember again that the 
aluminum panel is about 1 h times heavier than the composite 
panels. Therefore, utilizing composite panels can greatly 
increase the strength-to-weight ratios. 

The next composite study varies the lamination angle (0,- 
<p,<p) of a three-layer symmetric angle-ply laminate. Again, 
two cases were run to compare transmission loss data. The 
first case uses the same data as the simply supported case 
above. Figure 4.26 shows the composite panel data versus the 
brass panel at 80 Hz. For the square panel, note that all 
transmission loss data for the composite panel are higher than 
the data for the brass panel. The uncoupled natural frequency 
is close to the coupled frequency for the brass panel. Thus 
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the output pressure is higher than the input pressure and the 
transmission loss values are negative. Remember, the 
transmission loss characteristics are totally dependent on the 
input frequency. 

The second case uses the same data as the clamped case 
above. Figure 4.27 shows the comparison of the transmission 
loss results between the composite panel and the aluminum 
panel. Both panels have a dimensional ratio of 2:1 and the 
results show that in the case of the composite panel, the 
transmission loss is highly dependent on the lamination angle 
of the different layers. 

The results above for the composite panels are simple 
examples showing the coupling of the plate/cavity system. At 
this point, the design process becomes critical in selecting 
the correct composite laminate. Composites offer many design 
variables, such as lamination angle, layer thickness and 
volume-fraction of graphite fibers. The designer must know 
the operating frequency range of the vibrating structure 
before deciding which laminate to use. Using this coupled 
BEM-FEM technique, along with the design variables, an optimal 
composite laminate can be chosen, which will offer excellent 
transmission loss characteristics with a low weight advantage. 

4.3 Coupled Random Response of Structural- 
Acoustic Problems 

As shown in the harmonic response problems, modelling the 
coupling of a system becomes very important. For cavity- 
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Figure 4.26 Transmission Loss of Simply Supported, 
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Figure 4.27 Transmission Loss of Clamped, Three-Layer 
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backed plate problems, as shown in Figure 3.3, the first 
coupled natural frequency of the system is increased. The 
same problems from Section 4.2 are now analyzed with random 
white noise excitation. Isotropic and composite panels with 
the same structure and acoustic characteristics are modelled. 
The coupled random response is compared with the uncoupled 
structure response to show the effect of coupling. The 
coupled system involves analysis of the duct and the structure 
together whereas the uncoupled system neglects the effect of 
the cavity. 

First, the system with a simply supported panel is 
modelled. The square brass panel is 20 cm. by 20 cm. and 
0.09144 cm. thick. Figure 4.28 shows the random response of 
the brass panel for the coupled and the uncoupled systems. 
The input sound pressure level is varied from 90 to 150 dB. 
The maximum Root Mean Square (RMS) value of deflection (W,^) 
at 150 dB is approximately half of the thickness. As 
expected, the coupled system produces lower RMS (W^J values, 
since the first coupled natural frequency is increased. 
Similar to Equation (3.49), the response is dependent on the 
inverse of o) n 3 and thus the RMS response is increased since the 
first coupled natural frequency is increased. Figure 4.29 
models the same simply supported plate and cavity setup, but 
this time the panel is made of graphite/ epoxy . The duct 
cavity characteristics and brass panel properties are given on 
page 84 and the seven-layer Graphite/Epoxy panel 
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Figure 4.28 Random Response of Simply Supported Brass 
Panel and Acoustic Duct System. 
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Figure 4.29 Random Response of Simply Supported Composite 
Panel and Acoustic Duct System. 
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( 079070 o / 907079070 o ) properties are given on page 92. As in 
the earlier case, 25 finite elements and 150 boundary elements 
are used. Once again, the coupled system produces lower RMS 
( w m«) values. Note that the composite panel response is 
affected to a greater extent by the coupling than the brass 
panel. The composite panel is affected more, because the 
coupled natural frequency is increased from the uncoupled 
natural frequency by a greater amount. The uncoupled panel 
results agree with classical and other finite element 
solutions. 

The next coupled system consists of a clamped rectangular 
panel and the acoustic duct. The panel is 30.48 x 15.24 cm. 
and 0.16256 cm. thick. The duct cavity characteristics and 
aluminum panel properties are given on page 86 and the same 
Graphite/ Epoxy panel is used from above. Figure 4.30 shows 
the random response of the coupled and uncoupled clamped 
aluminum panel. Figure 4.31 shows the random response of a 
clamped composite panel. Once again, 25 finite elements and 
150 boundary elements are used. As expected, the RMS (W,,^) 
responses of the coupled systems are less than the responses 
of the uncoupled system. 

Figures 4.28 to 4.31 help us to infer that the coupling 
between sound and structure, greatly affects the random 
response of the system. 
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Figure 4.30 Random Response of Clamped Aluminum 
Panel and Acoustic Duct System. 
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Figure 4.31 Random Response of Clamped Composite 
Panel and Acoustic Duct System. 
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4.4 Modelling The Thermal Acoustic Fatigue Apparatus 
(TAFA) at NASA Langley Research Center 

Until now, we have seen that the coupled boundary/ finite 
element procedure gives accurate and reasonable results. 
Future high speed aircraft such as the National Aerospace 
Plane (NASP) , High Speed Civil Transport (HSCT) and the next 
generation of fighters will be subjected to intense thermal 
and acoustic environments [69], In an attempt to test these 
environments, the Thermal Acoustic Fatigue Apparatus (TAFA) at 
NASA Langley Research Center was built. TAFA, shown in Figure 
4.32, is a grazing incidence, high-intensity noise apparatus 
with capability of generating sound pressure levels of 168 dB 
and temperatures of 2000 °F. In this study, the temperature 
effects are not considered, but we account for acoustic 
pressures to be harmonic or random with a frequency input 
range of 0-500 Hz. 

4.4.1 Harmonic Response of Panels in TAFA Facility 

Three plates, brass, aluminum and graphite/epoxy, of 
15x15x0.036 in. are considered. The plate and the acoustic 
test section are discretized into 9 finite elements and 384 
boundary elements, respectively. The input and output 
pressure levels are specified as boundary conditions. The 
output boundary condition is set at 60 dB and the input 
pressure level varies from 90 to 140 dB. The input frequency 
of the external source is set at 20 Hz, because this is close 
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Clamping 





to the first natural frequency of the brass panel. This 
excites a first mode shape of each plate. The center 
displacements of the simply supported plates of the coupled 
system versus the input sound pressure level are shown in 
Figure 4.33. The composite panel is made of graphite/epoxy 
with a [ 4 5°/ -4 5°/ 4 5°/ -4 5°], layup. The composite panel has the 
following properties: E,,=22.5E6 psi, E 22 =1.17E6 psi, v n =.22, 
G 12 =0.66E6 psi, p=0.056 lb/ in 3 and ft =0.01. 

Figure 4.34 shows the displacement results for the 
clamped aluminum, brass and composite panels. From the 
Figures 4.33 and 4.34, we notice that the deflections of the 
composite panel are lower or equivalent to the deflections of 
the aluminum and brass panels. The considerable advantage of 
the composite panel is in the strength-to-weight ratio. The 
weights of the composite, aluminum and brass panels are .45 
lb, .81 lb and 2.49 lb, respectively. 

Two other features of the TAFA facility were also 
analyzed. First, the center panel deflection was calculated 
as a function of frequency for the simply supported composite 
panel listed above. Secondly, the pressure was analyzed 1/2" 
from the center of the panel on the interior of the duct. The 
center panel deflection and the acoustic pressures, as well as 
the uncoupled acoustic response, are seen in Figure 4.35. The 
acoustic response is defined as a sound pressure level, 


106 



Wmax/h 



— e- Brass Plate — ' *— Aluminum Plate Composite Plate 

Figure 4.33 Harmonic Response of Simply Supported Panels. 
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Figure 4.34 Harmonic Response of Clamped Panels. 
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Sound Pressure Level =20 Log 


(4.7) 


r*f 


where P ref =2xlO' 5 Pa. From Figure 4.35, we can see that the 
first coupled natural frequency of the panel is around 42 Hz. 
As the panel resonates at its natural frequencies, large 
displacements are occurring in the panel and the acoustic 
domain. These large acoustic displacements produce higher 
acoustic pressures and hence lead to the higher acoustic 
response peaks. 

In order to find the uncoupled acoustic natural 
frequencies, the same analysis was run in the TAFA facility 
with a rigid wall in place of the panel. The first three 
acoustic natural frequencies are 109 Hz, 179 Hz, and 240 Hz. 
The second uncoupled plate natural frequency was 98.7 Hz. 
Figure 4.36 shows the structure and acoustic response plotted 
at x=12 . 5 in. and y=7 . 5 in. Figure 4.3 6 shows that the 
structure is driving the system at the second natural 
frequency of the panel (98.7 Hz). Also from Figure 4.35, we 
can see that the acoustic response curve mirrors the structure 
response curve. This figure shows the strength of the coupled 
BEM-FEM formulation. The acoustic and structure response can 
simultaneously be monitored in the coupled system. 

The mode shape of the panel and the distribution of 
acoustic pressure are shown for the first coupled natural 
frequency (42 Hz) in Figure 4.37. Figure 4.38 shows the mode 
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Figure 4.35 Structure and Acoustic Response for 

Simply Supported Composite Panel in TAFA. 
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Figure 4.36 Structure and Acoustic Response for 

Simply Supported Composite Panel in TAFA 
at x=12.5" and y=7.5 M . 
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shape and the pressure distribution at 100 Hz, which is close 
to the second coupled natural frequency. Figure 4.38 shows 
that the panel is vibrating similar to a 2:1 mode shape. 

4.4.2 Random Response of Panels in TAFA Facility 

A random response of the coupled system is now analyzed, 
since the TAFA facility can produce a random acoustic input. 
A cut-off frequency of 400 Hz is assumed, which easily covers 
the first few modes of the composite panel. The same 
composite panel from Section 4.4.1 is used again. Figure 4.39 
shows the random response for the simply supported and the 
clamped composite panels, where the output noise level is set 
at 60 dB and the input noise level is varied from 130 to 160 
dB. The overall sound spectral level in Figure 4.39 is the 
difference between the input noise and the output noise. The 
random response of the clamped composite panel is compared 
with the response of the simply supported composite panel in 
Figure 4.39. The coupled RMS displacement of the clamped 
panel is less than the displacement of the simply supported 
panel . 
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Figure 4.37 Panel Mode Shape (Top) and Acoustic Pressure 
Distribution (Bottom) at 42 Hz. 
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Figure 4.38 Panel Mode Shape (Top) and Acoustic Pressure 
Distribution (Bottom) at 100 Hz. 
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Figure 4.39 Random Response of Clamped and 
Simply Supported Panels in the 
TAFA Facility. 





Chapter 5 


Summary and Conclusions 

5.1 Concluding Remarks 

In this study, we have presented the principles of a very 
useful tool for analyzing real world, sound-structure 
environments. In order to accurately model these complex 
coupled problems, coupling of the structure and acoustic 
environments must be done carefully. The coupled system 
effects were shown to be very important for certain types of 
problems. For example, the cavity-backed plate problem showed 
a significant first frequency shift due to the coupling. In 
our attempt to model real environments, random response 
analysis of the coupled system was also performed. The 
coupled harmonic response gave very accurate results when 
compared with existing analytical, exact and experimental data 
for simplistic problems. The Thermal Acoustic Fatigue 
Apparatus at NASA Langley Research Center was modelled and 
analyzed. The coupled BEM-FEM method shows the relationship 
between the interior acoustic pressure and the structure 
response for a 0-400 Hz frequency range. 

The boundary element method was also applied to interior 
acoustic domain problems. Two and three-dimensional analyses 
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were performed on simple and complex geometric-shaped ducts. 
Constant and linear boundary elements gave very satisfactory 
results for the 2-D case and constant elements for the 3-D 
case. 

Currently, industry is striving to couple the boundary 
element and finite element methods in order to more 
effectively model aerospace and automobile systems. This 
study is an attempt to enhance our understanding in that 
direction. Applying random analysis and adapting to composite 
materials allows designers to model real world problems. The 
random response formulation of the coupled finite element and 
boundary element method can provide a tremendous approximation 
method for complex coupled structural-acoustic systems. 

5.2 Future Work 

Many areas of research still need to be considered in 
order to more accurately model real world systems. In order 
to continue this study, four areas of research can be 
completed in the near future. First, the three-dimensional 
constant boundary element needs to be extended to a linear or 
even a quadratic type element. This will allow users to model 
more complex geometry systems as well as improve the accuracy 
of the results. Secondly, temperature effects can also be 
considered for the acoustic and the structural domains. 
Temperature effects will greatly affect the structure, which 
will in— turn greatly affect the acoustic response. The next 


116 



area of research could be in the harmonic nonlinear area. The 
structure can be modelled with nonlinear analysis where the 
time components of the nonlinear terms neglect the period 
doubling effects. The last area of interest is in the random 
nonlinear analysis. The nonlinear analysis of the structure 
can be incorporated into the current finite element analysis 
of the structure. 
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APPENDIX A 


The integrals in Equations (2.18), (2.22) and (2.25) must 
be evaluated in order to obtain each element of the influence 
matrices, G mj and The easiest and most efficient way of 

solving these integrals is approximating them with the Gauss 
quadrature scheme. The first step using Gauss quadrature is 
to approximate the integral by transferring the limits from -1 
to 1 and then apply the Gauss rule to the new integral. In 
other words, 

a l nn 

I=\f(x)dx = f F(xJdx t =i=£w M ,F(xJ , (A.l) 

i 'i "-i 

where W w are the weighting functions, x w are the gauss points, 
and nn is the total number of gauss points. A four-point 
gauss scheme is used and the weight and point values are shown 
below, 

Table A.l Gauss Points and Gauss Weights 


Point 

Gauss Point 

Weight 




W w 


1,2 

±.86113631 

.34785485 


3,4 

±.33998104 

.65214515 



Note that the transfer from the exact integral to the 
Gauss integral requires a transformation in the coordinates 
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and hence a scaling factor in terras of the element size is 
produced. A four-point Gauss scheme was found to be 
sufficient in calculating the integrals of Equations (2.18), 
(2.22) and (2.25) . 

The Struve functions given in Equations (2.35) and (2.36) 
are not defined as Fortran functions or in the IMSL libraries. 
Therefore, approximations to the Struve functions must be 
incorporated. The approximation to the Struve function of the 
first kind is 



^2k*\ ( Z ) 
2/C + 1 ' 


(A. 2) 


or 


( % ) ~ — 


n [ 1 


z) i7 3 ( Z ) Jj (z) J 7 (z) 


(A. 3) 


where Jj is the first Bessel function of the ith kind. The 
approximation to the Struve function of the second kind is 


S T i ( z ) ~ — — ‘J 0 ( z ) + — zZ 


tt n 


7 T ^ 


Jk(z) 

4k 2 -l 


(A. 4) 


or 


77 n n oK n ] 3 15 35 


)- 


(A. 5) 


These approximations to the Struve functions give very 
accurate solutions [66]. 
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APPENDIX B 


The element selected for this study is a four-node 
rectangular element as shown in Figure 3.1 [67]. The element 
is a C 1 conforming element with twenty-four degrees of freedom 
(six at each node w, w x/ w y , w^, u, and v) . The shape 
functions for the transverse and in-plane displacements are 

H w ={ 1 x y x 2 xy y 2 x 3 x 2 y xy 2 y 3 

(B. i; 

x 3 y x 2 y 2 xy 3 x 3 y 2 x 2 y 3 x 3 y 3 } , 

and 

H u ={ 1 x y xy 0 0 0 0}, 

H v ={ 0 0 0 0 1 x y xy} . (B.2) 

The generalized coordinates from Equations (3.4), (3.5) and 

(3.6) are defined as 

(a, a 2 0:3 a 4 a 5 a 6 a 7 a 8 

a = 

a 9 a i0 a il a i2 a !3 a i4 a i5 a i6^ / 

and 

0 = {/3, 0 7 j 8 3 jS 4 )8 5 /3 6 /3 7 /3 g > . (B.3) 
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The degrees of freedom in bending and membrane are 


i W b} T = i W t W 2 W 3 V A W xl' * W x4 W >r ' W y4 W xyl' ' . ( 

and 

u 2 u, u 4 V, VjV,V 4 ]. ( 


The inverse of the transformation matrix, [T b ] 
Equation (3.6) is given as 


1 

0 

0 

0 

0 

0 

0 

0 

1 

al 

0 

al 2 
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0 

ai 3 

0 

1 

al 

bl 

al 2 

al bl 

bl 7 

al 3 

ai 2 bl 

1 

0 

bl 

0 

0 

bl 2 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

1 

0 

2al 

0 

0 

3 al 2 

0 

0 

1 

0 

2 al 

bl 

0 

3al 2 

2 ai bl 

0 

1 

0 

0 

bl 

0 

0 

0 

0 

0 

l 

0 

0 

0 

0 

0 

0 
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al 
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0 

0 
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0 

0 
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B.4) 


B. 5) 


from 
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(B.6) 

The transformation matrix [T m ] is given by 


[T m ] 


1 

-1 

al 

l 


o 
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al 

-1 
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al bl al bl al bl al bl 

(B.7) 

where al and bl are the length and width of the element. 
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From Equation (3.8), the matrix [C b ] is given as 


[C b ] = 


-if 



0002006* 2 y 
000002 0 0 
000020 0 4x 


0 

0 

6 xy 

2y 2 

0 

6 xy 2 

2y 3 

6xy 3 “ 


2x 

e y 

0 

2X 2 

6 xy 

2* 3 

6x 2 y 

6x 3 y 

(B.8) 

4y 

0 

6x 2 

8 xy 

6y 2 

12x 2 y 

12xy 2 

18x 2 y 2 



The matrix [C m ] from Equation (3.9) is given as 



H 


- 

[Cm ] = 

u, a 

- 

OlOyOOOO 

OOOOOOlx 




OOlxOlOy 


The stress-strain relationship in Equation (3.10) 
involves the reduced stiffness matrix, [Q] . The Q, 6 , and Q 26 
terms are set to zero, due to orthogonality. Each element of 
the [Q] matrix is a function of the composite material 
characteristics including: Young's modulus (E u and E 22 ) , 

Poisson's ratio {v n and v 2x ) , and the shear modulus (G 12 ) . The 
[Q] matrix components are defined as: 
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Q.r 




Qn~ 


-22 




0)2 


y i 2 ~^22 
1 -^ 12^21 


C?66 = G,2- ( B * 10 ) 

If a layer is rotated by an angle 0, then the reduced 
stiffness matrix must be transformed to the [Q] reduced 
stiffness matrix. The transformation matrix, [T] , is given by 


[T] 


cos 2 0 sin 2 0 2sin0 cos0 

sin 2 0 cos 2 0 -2sin0 cos0 

l-sin0 cos0 sin0 cos0 cos 2 0-sin 2 0 


(B. 11) 


The transformed reduced stiffness matrix, [Q] , is calculated 
as 

[0] = [T]-‘ [ Q ] [T]- r . ( B . 12 ) 

The terms of the transformed stiffness matrix are given as 


On = On cos 4 0 + 2 (,Q n +2Q 66 ) sin 2 0 cos 2 0+O 22 sin 4 0 


Q ,2 = ( 0 n + C ) 22 - 4 C ) 66 ) sin 2 0 cos 2 0 +O 12 (sin 4 0 +cos 4 0 ) 
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0 22 = O,,sin 4 0 + 2 (0,2+20^) sin 2 0 cos 2 0 +£> 22 cos 4 0 


0,6 = (011-012-2066)51110 cos 3 0+(O,2-O 22 +2O66) sin 3 0 cos 0 


0 2 6 = (O„-O, 2 -2O66)sin 3 0 cos0+ (0 I2 -0 22 + 2Q w ) sin0 cos 3 0 


0e6 = (0 u+0 22 -20i 2 -2Q 66 ) sin 2 0 cos 2 0+Q w (sin 4 0+cos 4 0) • 


(B . 13 ) 
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